Functional renormalization group and variational Monte Carlo studies of the 
electronic instabilities in graphene near 1/4 doping 
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We study the electronic instabilities of near 1 /4 electron doped graphene using the singular- mode 
functional renormalization group, with a self-adaptive k-mesh to improve the treatment of the van 
Hove singularities, and variational Monte-Carlo method. At 1/4 doping the system is a chiral spin 
density wave state exhibiting the anomalous quantized Hall effect. When the doping deviates from 
1/4, the d^2_y2 +idxy Cooper pairing becomes the leading instability. Our results suggest that near 
1/4 electron- or hole- doping (away from the neutral point) the graphene is either a Chern insulator 
or a topoligical superconductor. 

PACS numbers: 74.20.-z, 74.20.Rp, 74.70.Wz, 81.05.ue, 71.27.-|-a 



INTRODUCTION 

Graphene, a single atomic layer of graphite, has been 
a focus of interest since the pioneering work of Novoselov 
and Gcim[l]. At the fundamental level the past research 
activities on graphene mostly focused on exploring the 
consequences of the unique Dirac-like bandstructure[2]. 
On the experimental side, few exceptions include the ob- 
servation of the fractional quantum Hall effect [3, 4], the 
detection of the Fermi velocity renormalization[5], and 
the possible observation of "plasmaron" in angle-rcsovled 
photosemission[6]. In general the effects of electron- 
electron interaction on the properties of graphene re- 
main a frontier of this field. Previously based on the 
resonating- valence-bond [7] concept Pathak et al.[8] and 
Black-Schaffer and Doniach [9] proposed that doped 
graphene should be a high temperature superconduc- 
tor with d + id' pairing symmetry. (Henceforth d and 
d' are used to denote interchangeably d^^-y^ and d^y 
symmetries, respectively.) In particular, the possibility 
of unusual superconductivity and other orders in doped 
graphene with van Hove singularities at (or near) the 
Fermi level becomes a hot issue. [10, 11] Most recently by 
a perturbative renormalization group calculation Nandk- 
ishore et al. concluded that the van Hove singularities on 
the Fermi surface drive chiral d + id' superconductivity 
in the limit of vanishing interaction strength[12]. 

On a different front Tao Li recently proposed that due 
to the existence of Fermi surface nesting 1/4 electron 
doped Hubbard model on honeycomb lattice favors the 
formation of a magnetic insulating state which possesses 
nonzero spin chirality and exhibit the anomalous quan- 
tized Hall effect, hence is a Chern insulator[13]. Thus 
near quarter doping graphene suddenly becomes a play- 
ing ground where either a Chern insulator or a topolog- 
ical superconductor can potentially be realized. Because 
the realization of either phases in heavily doped graphene 



will be truly exciting, we feel it is meaningful to examine 
this problem using the more realistic band structure and 
interaction parameters. 

In view of the heavy doping we use the Hubbard inter- 
action to model the screened Coulomb interaction. We 
perform singular-mode functional renormalization group 
(SM-FRG)[14] and variational Monte Carlo (VMC) cal- 
culations to address the possible electronic instabilities. 
Since the interaction strength is estimated to be a frac- 
tion of the band width, we believe SM-FRG should yield 
qualitatively correct answer. The VMC is used to further 
confirm such belief. The main results are summarized 
as follows. At 1/4 electron doping and for interaction 
strength appropriate for graphene we found the chiral 
spin density wave (SDW) state is the dominating insta- 
bility. When the doping level slightly deviates from 1/4 
we find the d + id' pairing instability surpasses that of 
the chiral SDW. We propose a schematic phase diagram 
in Fig. 6(b). As in pnictides and ovcrdoped cupratcs[15], 
the pairing mechanism is due to a strong scattering chan- 
nel shared by the SDW and pairing. 



MODEL 

The real-space hamiltonian we used is given by 

(ij)a I 

+ ]^yYn,n^+S, (1) 

iS 

where (ij) denotes bonds connecting sites i and j, a is 
the spin polarity, fj, is the chemical potential, A'e is the 
total electron number operator, the [/-term is the on-site 
Hubbard interaction and V is the Coulomb interaction 
on nearest-neighbor bonds 5. The honeycomb lattice has 
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two sublattices which we denote as A and B henceforth. 
As suggested in Ref. [2] we take ti = 2.8 eV, t2 =0.1 eV, 
and ts = 0.07 eV for hoppings between the first, second 
the third neighbors, respectively, and set U = 3.6<i. As 
for V, we expect V < U in doped graphene, and take 
V ^ ti as a typical upper bound. Theoretically, enriched 
phases may appear for even larger values oi V/U.[ll, 16] 



METHOD 

The SM-FRG niethod[14] wc used is a modification 
of the FRG method[17] applied to the cuprates[18] and 
pnictides[19]. Fig.l (a) shows a generic 4-point ver- 
tex function ri234 which appears in the interaction 
c\cl{—T 1234) C3C4. Here 1, 2, 3, 4 represent momentum (or 
real space position) and sublatticc. The spin a and r arc 
conserved during fcrmion propagation, and will be sup- 
pressed henceforth. Figs.l(b)-(d) are rearrangements of 
(a) into the pairing (P), the crossing (C) and the direct 
(D) channels in such a way that a collective momentum q 
can be associated and the other momentum dependence 
can be decomposed as, 

i k+q,-k,-p,p+q 

mn 

^ k+q,p,k,p+q 

mn 

^ k+q,p,p+q.k 

mn 

Here {/m} is a set of orthonormal lattice form fac- 
tors. For honeycomb lattice the form factor label m 
also includes a sublattice label, m = (m, a) with a — 
A/B, within our choice of C^y point group with re- 
spect to the atomic site. (See Appendix.) The de- 
composition into each channel would be exact if the 
set is complete. In practice, however, a set of a few 
form factors is often sufficient to capture the symme- 
try of the order parameters associated with leading 
instabilities. [14] The momentum space form factors arc 
the Fourier transform of the real-space ones: 1) on-site, 
fso{r) = 1; 2) 1st neighbor bonds, /^^(r) = 
fddr) = y^cos{W,) and fd,,{r) = y/2/3sm{W,), 
where / = 2 and 9^ is the azimuthal angle of r;[20] 3) 
2nd neighbor bonds, /s2(r) = ^^1/6, fp^ = \/l/3cos 0r, 
fp^,{r) = yiTSsin^r, /d.(r) = v^cos2^r, /d„ (r) = 
■\/l/3sin20r, ff^{v) = a/1/6cos Wr. These form factors 
(combined with sublattice labels) are used in all channels 
on equal footing. We have tested that further neighbor 
form factors do not change the results qualitatively. 

The one-loop correction to the flow of the vertex func- 
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FIG. 1: A generic 4-point vertex (a) is rearranged into the 
pairing (b), crossing (c) and direct (d) channels. Here k, q, p 
are momenta, a and r denote spins which are conserved dur- 
ing fermion propagation, and m, n denote the form factor (see 
the text for details). The open arrows indicate collective prop- 
agation. 

tion can be written as, in matrix form, 

dP/dK = Px;p^, 

dc/dk = cx;„c, 

dD/dK^{C~D)^^^D + Dx'pHiC-D), (3) 

where the collective momentum q is left implicit for 
brevity, x'pp/ph ioo^i integrations projected by the 
form factors (See Appendix for details), and A is the 
running cutoff energy. Integrating Eq. (3) with respect 
to A yields the ladder approximation. However, since 
dP. dC and dD add up to the full change dF, the full 
flow equations for P, C and D should be given by 

dP/dk = dP/dk + P{dC/dk + dD/dk), 
dC/dk = dC/dk + C{dP/dk + dD/dk), 
dD/dk = dD/dk + D{dP/dk + dC/dk), (4) 

where the P, C and D are the projection operators in 
the sense of Eq. (2), and we used the fact that K[dK) = 
OK for K = P,C, D. In Eq. (4) the terms preceded 
by the projection operators represent the overlaps of the 
three different channels. It is those terms which allow 
pairing to be induced by virtual particle-hole scattering 
processes[15]. 

It can be shown that the effective interaction in the 
superconducting (SC), spin density wave (SDW) and 
charge density wave (CDW) channels are given by Vsc = 
-P, Vsdw = C, and Vcdw = C - 2D, respectively. By 
singular value decomposition, we determine the leading 
instability in each channel, 

Vr{<lx) = E S^rximH^in), (5) 

a 

where X = sc, sdw, cdw, Sx is the singular value of the 
a-th singular mode, (jj^ and i/'x ^^'^ right and left 
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FIG. 2: (Color online) Results for 5 = 1/4 and V = 0. (a) 
The Fermi surface. The brightness on the surface represents 
the momentum space gap function associated with one of the 
degenerate pairing modes, (b) FRG flow of (the inverse of) 
the most negative singular values S in the SC (blue solid 
line), SDW (green dot-dashed line) and CDW (red dashed 
line) channels, (c) and (d) are the renormalized interaction 
K"?™ for = {so, A), and 1/™™ for m = (di, A), as functions 
of the collective momentum q. The hexagons in (a), (c) and 
(d) indicate the Brillouine zone boundary. 



eigen vectors of Vx , respectively. We fix the phase of the 
eigen vectors by requiring Re[X]m 4'x(''^)^xi^''-)] > so 
that Sx < corresponds to an attractive mode in the 
X-channel. In the pairing channel qsc = addresses 
the Cooper instability. The pairing function in the sub- 
lattice basis can be constructed from 0"^, and a fur- 
ther unitary transform is needed to get the pairing func- 
tion in the band basis. (Sec the Appendix for details.) 
The ordering wave vector in the SDW/CDW channel 
q = Isdw/cdw is chosen at which \Vsdw/cdwi^)\ is niax- 
imal. We note that such a vector has symmetry-related 
images, and may change during the FRG How before set- 
tling down to fixed values. The RG flow is stopped if any 
of \P\,nax, \C\rnax, Or \D\„iax becomcs roughly 10 times 
of the bandwidth. [21] 

More technical details can be found in the Appendix. 



RESULTS AND DISCUSSION 

We define the doping level hy 5 = rie — 1 where rie 
is the number of electrons per site. Wc first discuss the 
results for (5=1/4 and V ^0. Fig. 2(a) shows the Fermi 
surface in this case, which is well nested and close to 
the van Hove singularities (the mid points of edges of 



the outer hexagon). The flow of the most negative sin- 
gular values (denoted as S) in the SC, SDW and CDW 
channels are shown in Fig. 2(b). Clearly the SDW (green 
dot-dashed line) is the leading instability. This is be- 
cause the SDW scattering is already attractive at high 
energies and is further enhanced by the Fermi surface 
nesting shown in Fig. 2(a) down to the lowest energies. 
This SDW singular mode 4>sdw {rn) has a dominant value 
for m — {so, A/ B), showing that the magnetic order- 
ing moment is site-local. It is also identifiable in the 
renormalized interaction ^^^^(q) for m = {so,A) shown 
in Fig. 2(c), which has strong peaks at three indepen- 
dent momenta Qi = (0, 27r/-\/3), Q2 = (— tt, 7r/-\/3), 
Q3 = (7r,7r/-\/3) and their symmetric images. They 
define the possible ordering vectors qsdw for the SDW 
order. In contrast, the scattering associated with pair- 
ing (blue solid line) is initially repulsive in all channels, 
and only becomes attractive in the d-wave channel af- 
ter the SDW scattering grows strong. Regarding Cooper 
pairing we find two degenerate leading form factors: the 
dx'^-y'^ and d^y doublets. One of these is used to gener- 
ate the momentum space gap function shown in Fig. 2(a). 
The singular pairing mode (j)sc{fn) is nonzero for the 2nd 
neighbor bonds, but the amplitude is about one order of 
magnitude smaller than that for the 1st neighbor bonds, 
hence justifying the cutoff in the real-space range of the 
form factors. The renormalized Vy^'"(q) for m = {di,A) 
in Fig. 2(d) has negative (but weak) peaks at q = 0, con- 
firming the Cooper instability at this momentum. The 
CDW channel (red dashed line) also becomes weakly at- 
tractive from Fig. 2(b). We checked that the singular 
mode 4>cdw{fn) has significant values for m — {si,A/B) 
and m = (di^i', A/i?), showing that it is an extended 
CDW. The mixture of si and di^y is due to the fact that 
the CDW wave vector q = (27r/3, tt) (or its symmetry im- 
ages) is not invariant under the point group operations. 
However, the CDW channel remains weak in our case 
hence will not be discussed in the rest of the paper. (The 
merging of SC and CDW lines in Fig. 2(b) is induced by 
the diverging SDW channel via the overlapping among 
these channels.) 

The above results indicate three independent and de- 
generate SDW momenta (apart from the global spin 
SU{2) symmetry). A calculation that keeps the 
symmetry-breaking self-energy flow is needed to de- 
cide which combination of them is realized in the or- 
dered state, but this is beyond the scope of the present 
work. Alternatively, one may resort to a Landau the- 
ory or mean field theory. Indeed, according to the mean 
field theory of Ref. [13], a particular linear combination 
[13] (Sr,a) = Mae*"^^-^ + Mie'^^ '^ + Mae**^^'^ and 
(SR,i3) = Mae**^^-^ - Mie'^Ji-^ - Mae^^^^-^, gives the 
most energetically favorable spin structure. Here R la- 
bels unit cell and Mi_2,3 are three mutually orthogonal 
and equal length vectors. The handedness of the Mi. 2, 3 
triad breaks time-reversal and spatial reficction symmc- 
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FIG. 3: (Color online) Chiral SDW c 
comb lattice, (a) The spins on the blac 
sublattices (of different gray scales) ai 
-Ml - M2 + Ma, Ml - M2 - M3, -M 
tively. (b) A 3D perspective view of the 



try. The resulting four-sublattice cl 
shown in Fig. 3. 

Since our result is at odd with that of Rcf. [12] which 
applies in the limit of vanishing interaction strength, 
we further check the above conclusion by a variational 
Monte-Carlo (VMC) calculation using exactly the same 
parameters as for Figs. 2. We adopted the partially- 
projected mean-field wave-functions[22] as our trial wave- 
functions, with the the form factors guided by the present 
SM-FRG result. Fig. 4 shows the energy gain per site 
due to d 4- id! pairing on 12 x 12 (circles) and 18 x 18 
(triangles) lattices, showing negligible size dependence. 
Wc then compare to the energy gain associated with the 
chiral SDW (squares). It is clear that the SDW state is 
far more energetically favorable than the chiral d + id'- 
SC state at this doping level. The reason that the SDW 
state is realized in our lattice model lies in the fact that 
the perfect Fermi surface nesting is as important as the 
inter-saddle scattering addressed in Ref.[12]. 

Below 1/4 doping, the the Fermi level moves away from 
the van Hove points and the Fermi surface nesting wors- 
ens. This is shown in Fig. 5(a) for 5 = 0.211 as an ex- 
ample. The bare interactions are still set as U = 3.6ii 
and V ^ Q. The SDW scattering is still attractive at 
high energies. As seen in Eq.(4), this relatively strong 
SDW channel causes attraction in the SC channel via 
overlap between these channels (terms with the projec- 
tion operator in Eq. (4)). At even lower energy scales 
the pairing channel attraction (with q = 0) continues to 
grow due to the Cooper instability, while the enhance- 
ment of the SDW scattering is saturated due to the lack 
of Fermi surface nesting. As the result the pairing insta- 
bility surpasses the SDW instability at the lowest energy 
scale. This is shown by the flow of the singular values 
in Fig. 5(b). It is worth to mention that precisely the 
same phenomenon was observed in the FRG studies of 
the cuprates and pnictides[15, 18]. A close inspection 
of the eigenvectors (j)sc{m) associated with the most di- 
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FIG. 4: (Color online) Variational Monte Carlo results for 
the energy gain per site, AiJ, versus the variational order 
parameters A for the d + id'-SC (circles and triangles) and 
the chiral SDW states (squares) at the doping level (5 = 1/4. 
The lattice sizes are given in the legends. Lines are guides to 
the eye. 
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FIG. 5: (Color online) The same plots as in Fig. 2 but for 
5 = 0.211. Note the splitting of the SDW peaks in panel (c) 
signifies the incommensurate SDW instability. 

verging superconducting pairing channel again find the 
degenerate d^2_y2 and dxy doublets, with dominant am- 
plitudes for m ~ (di^i' , A/ B). The momentum space 
gap function of one of the pairing modes is shown in 
Fig. 5(a). Fig. 5(c) shows the renormalized V^^{q) for 
TO = {so,A), which shows weak peaks at six indepen- 
dent and incommensurate momenta. Fig. 5(d) shows the 
renormalized Vj,™™(q) for to = {di,A), which shows a 
strong negative peak at q = 0. 

The above results imply degenerate c?-wavc pairing 



5 




0.15 0.2 0.25 0.3 1/4 

5 8 



FIG. 6: (Color online) (a) The FRG diverging energy scale 
Ac plotted as a function of doping level near 5 = 1/4. crosses 
and open circles represent Ac associated with the SC and 
SDW channel, respectively. V = 0,ti for solid and dashed 
lines, respectively, (b) A schematic temperature-doping phase 
diagram near 5 = 1/4 in linear scales. The grey region denotes 
the transition between SC and SDW. 

instabilities. As for the degenerate SDW instabilities, 
additional analysis, such as the mean field theory or 
Ginzburg-Landau theory, is needed to fix the structure of 
the pairing function in the ordered state. To a large ex- 
tent, this kind of analysis has been performed in Ref. [12]. 
We have also performed simple mean field calculations 
using the renormalized pairing interaction. The result is 
that a time-reversal breaking dx2_y2 ± idxy-^&'ve pairing 
is always more favorable. This could have been antici- 
pated since both dx2-y2 and d^y form factors have nodes 
on the Fermi surface, a natural way to gain energy is to 
form the above chiral d-wave pairing, which gaps out the 
entire Fermi surface. 

We have also checked nearby doping and analyzed the 
competition between the incommensurate SDW and the 
SC state. In Fig. 6(a) we plot the higher diverging scale 
among these two competing orders as a function of dop- 
ing (solid line). We found a similar phase diagram near 
-1/4 (hole) doping (not shown), mirroring that of elec- 
tron doping. (Notice that the particle-hole symmetry is 
not exact in the presence of hopping integrals ^2.3)- 

We now discuss briefly the effect of the nearest neigh- 
bor interaction V. As a typical example, we set U = 3.6ii 
and V ^ ti, and perform the FRG calculations. We 
find that the results are qualitatively similar to the cases 
with V = 0, except that in the leading pairing singu- 
lar mode, (f>sc{d2,2' , A/ B) becomes slightly stronger, but 
still smaller than (l)sc{di,i' , A/ B) by a factor of 4 ^ 6. 



The phase diagram for V = ti is also drawn in Fig. 6(a) 
(dashed line). The critical scale is slightly higher than the 
case oi V = 0. In the SC region this is due to the slight 
enhancement of d-wave pairing on 2nd neighbor bonds. 
Unlike that claimed in Ref. [23], in all cases studied in 
this paper the /-wave pairing is not a leading instability. 

We end by presenting Fig. 6(b) as a schematic phase 
diagram in the temperature-doping plane. In reality 
when the doping level slightly deviates from 1/4, the ex- 
tra charges will be localized by the presence of disorder, 
which enables the system to stay in the Chern insulator 
state for a finite doping interval. In the transition region 
marked by gray, where the doped charges delocalize, in- 
commensurate SDW states with wave vectors shown in 
Fig. 5(c) will emerge. 



SUMMARY 

In summary, we have performed SM-FRG calculations 
for parameters suitable for graphene. Our results indi- 
cate that near 1/4 electron- or hole- doping, graphene is 
either a Chern insulator or a chiral d-wave superconduc- 
tors. Both phases are topological in nature, and deserve 
experimental efforts in searching them. 
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APPENDIX 

To illustrate the idea of the method, we first ignore the 
sublattice index, and return to it at a later stage. 

The one-loop contributions to the fiow of the irre- 
ducible 4-point vertex function is shown in Fig. 7, where 
(a) and (b) lead to the fiow of P and C, respectively, 
and (c)-(c) lead to the flow of D. The internal Greens 
functions arc convoluted with the form factors, hence 
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FIG. 7: One loop diagrams contributing to the flow of the 
the 4-point vertex function in the pairing channel (a), cross- 
ing channel (b), and direct channel (c)-(e). Here m,m'n,n' 
denote form factors, while the momentum and spin indices 
are left implicit. The open arrows indicate the flow of the col- 
lective momentum. The slashed lines are single-scale fermion 
propagators. The slash can be placed on either internal lines 
associated with the loop. 



where G is the bare fermion propagator, Sbz is the area 
of the Brillouine zone. Here A > is the infrared cut- 
off of the Matsubara frequency uJn ■ [24] As in usual FRG 
implementation, the self energy correction and frequency 
dependence of the vertex function are ignored. 

In general, the form factor /m(k) = J2r fm{^) exp(— ik- 
r), where fm(j) transforms according to an irreducible 
representation of the point group, and r is the relative 
position vector between the two fermion fields on each 
side of the diagrams in Fig.l (b)-(c). For two types of di- 
agrams to overlap, all of the four fermion fields sit within 
the range set by the form factors. Hence the projections 
in Eq.(4) are all preformed in real space. 

We now return to the honeycomb lattice with two 
sublattices. The necessary modifications are as follows: 
1) The sublattice label can be absorbed into the labels 
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FIG. 8: (Color online) Examples of self-adaptive k-mesh (a) 
and q-mesh (b) used for loop integrations and collective mo- 
menta for interactions, respectively. Notice that in (a) the 
mesh is too dense near the Fermi surface to be differentiated 
by naked eyes. 



1,2,3,4 in Figs.l, so in principle the label m in the 
form factor /m(r) also includes sublattice indices. How- 
ever, once r is fixed they are not independent. We ab- 
sorb an independent index into the form factor label, 
771 — > (m,a), where a labels fermion field 1 or 4 in 
Fig. 1(b), 1 or 4 in (c), and 1 or 3 in (d). This is rea- 
sonable because the point group operations do not mix 
the sublattice when the center is chosen to be atomic 
sites. 2) In the presence of sublattice indices the Green's 
functions arc matrices. 3) In order to ensure that in mo- 
mentum space P, G and D transform exactly as product 
of form factors, care must be taken in choosing the phase 
of the Bloch states for complex unit cells. 

The pairing function is determined as follows. A 
singular mode t/f^^ corresponds to a pairing operator 
Em=(m,a)clt(k)Cc('^)/m(k)*4^^^(-k) in the momen- 
tum space, where a™ is determined by 777 = (777, a) (for 
all associated vectors r) . The parity of the pairing matrix 
function under space inversion determines automatically 
whether it is a spin singlet or spin triplet. A further uni- 
tary transform can be used to get the momentum space 
gap function. 

In the current implementation of SM-FRG the sam- 
pling of momentum space (k and q) is performed on self- 
adaptive meshes. As an example. Fig. 8(a) shows the k- 
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mesh, which is progressively denser in approaching the 
Fermi surface or van Hove points (if they are close to the 
Fermi level). This is important for our problem because 
of the rapid variation of the Fermi velocity near the van 
Hove singularities. The k-mesh is obtained as follows. 
First define an energy scale V, (of the order of the band- 
width), and begin with 6 equal-area triangles spanning 
the BZ. Break a specific triangle into four smaller equal- 
area ones if any eigen energy |ek| < in the original tri- 
angle. Then lower the energy scale as f2 — ?> r2/& (6 > 1) 
and repeat the above process recursively. The centers of 
the triangles form the mesh points. Combining the trian- 
gle areas, they are used in the loop integrations in Eq.(6). 
In our implementation, we perform the above processes 
8-10 times, so that the last generation of triangles has a 
linear size of order 2~^tt ^ 2~^°7r, and the center of such 
triangles are sufficiently close to (or accidentally on) the 
Fermi surface. The q-mesh as in Fig. 8(b) used for the 
interactions includes all important scattering momenta: 
the origin and the high symmetry nesting vectors. We de- 
vise a function yyq such that it is zero at those important 
scattering momenta, and generate the mesh in a similar 
fashion as for k, except that 51 becomes an artificial scale 
and r]q is used in place of ek- Similar q-mcsh already 
appears in Rcf. [14]. 
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